<<<<<<< HEAD <<<<<<< HEAD ======= >>>>>>> d546939e94d695c6e9c44c085bc5ec2e2e14440e ======= >>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4
Ctrl+K

Lab Exercise 1

Lab Exercise 1#

The purpose of the first set of exercises is, on one hand, to familiarize with the programming environment of Python and Matlab, and on the other hand, to introduce the methods of representation and processing of telecommunication signals in these specific programming languages.

Part 1: Training#

Live Code

Press the following button to make python code interactive. It will connect you to a kernel once it says “ready” (might take a bit, especially the first time it runs).

Importing packages we will need later in Python#

from scipy import signal
# Ανατρέξτε στην τεκμηρίωση της βιβλιοθήκης scipy.signal
# https://docs.scipy.org/doc/scipy/reference/signal.html
from scipy.fft import fft, fftfreq
from scipy.fftpack import fftshift, ifftshift
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from numpy import random
import warnings
warnings.filterwarnings('ignore')
from ipywidgets import IntRangeSlider, widgets, Layout, HBox, VBox
from IPython.display import display, clear_output
print("Libraries added successfully!")

Create a scalar (one-dimensional) quantity#

s=2
print('s =',s)

s = 2;
disp(['s = ', num2str(s)]);

Create a vector of real values#

v=np.array([1,5,9])
print('v =',v)

v = [1,5,9];
disp(['v = ', mat2str(v)]);

Create a matrix of real values#

a=np.array([[1,2,3],[4,5,6],[7,8,9]])
print('a =',a)

a = [1,2,3; 4,5,6; 7,8,9];
disp(['a = ', mat2str(a)]);

Sum#

a+5

a_plus_5 = a + 5;
disp(a_plus_5);

Multiply#

b=s*v*2
print('b=',b)

b = s * v * 2;
disp(['b = ', mat2str(b)]);

Multiply element-wise#

np.multiply(v,b)

elementwise_product = v .* b;
disp(elementwise_product);

Check the length of a vector#

len(v)

length_v = length(v);
disp(length_v);

Check the size of a matrix#

a.shape   # για array: np.array(a.shape)

shape_a = size(a);  % This will return a 2-element vector: [rows, columns]
disp(shape_a);

Access specific elements of a matrix (Part 1)#

a[0,1]   # Η δεικτοδότηση αρχίζει από το 0

element_0_1 = a(1,2);  % MATLAB indexing starts at 1
disp(element_0_1);

Access specific elements of a matrix (Part 2)#

a[1,-1]   # Αρνητικές τιμές μετρούν από το τέλος, π.χ. το -1 αναφέρεται στο τελευταίο στοιχείο

element_1_end = a(2,end);  % 'end' keyword is used for the last element in MATLAB
disp(element_1_end);

Access a specific segment of a vector#

v1 = v[1:3]
v2 = v[1:2]
print('v1 =',v1)
print('v2 =',v2)   # ΠΡΟΣΟΧΗ: τα στοιχεία [2ο,3ο] δίνονται ως 1:3 και όχι ως 1:2

v1 = v(2:3);
v2 = v(2);
disp(['v1 = ', mat2str(v1)]);
disp(['v2 = ', mat2str(v2)]);

Access specific segments of a matrix#

a[0:2,:]   # Ομοίως: οι γραμμές 1 & 2 δίνονται ως 0:2 και όχι ως 0:1

a_slice = a(1:2,:);
disp(a_slice);

Create a vector with elements from 0 to 0.5 with a step of 0.1#

t=np.arange(0,0.5,0.1)
print('t=',t)

t = 0:0.1:0.5;  % The endpoint is inclusive in MATLAB, so you might want to adjust it if you want the exact same output as NumPy
t(end) = [];  % Remove the last element to match Python's exclusive endpoint
disp(['t = ', mat2str(t)]);

Part 2: Theory#

Sampling - Digitization#

Primary signals are mainly analog (continuous time). To represent and process them on our computer (or another digital machine), we must first digitize them. Consider a continuous-time signal \(x(t)\) with Fourier transform (Continuous Time Fourier Transform – CTFT):

\[ X(f)=\int_{-\infty}^{\infty} x(t)e^{-j2\pi ft} dt \]
# Parameters
fs = 1000  # Sampling frequency in Hz
T = 1 / fs  # Sampling period
L = 1000  # Number of samples
t = np.arange(0, L) * T  # Create time vector


# Output widget for the plots
graph_output0 = widgets.Output()

def update_plot(selected_frequencies_first):
    with graph_output0:
        clear_output(wait=True)  # Clear the previous plot
        # Generate the signal with the selected frequencies
        signal = np.zeros(len(t))
        for f in range(selected_frequencies_first[0], selected_frequencies_first[1] + 1):
            signal += np.sin(2 * np.pi * f * t)

        # Calculate the FFT of the original signal
        original_signal_fft = fft(signal)
        # Calculate frequencies for the FFT of the original signal
        original_freqs = fftfreq(L, T)[:L // 2]
        # Calculate magnitude of Fourier coefficients (amplitude) for the original signal
        original_magnitude = np.abs(original_signal_fft)[:L // 2]
        
        # Plotting
        fig, axs = plt.subplots(1, 2, figsize=(15, 5))
        
        # Original signal plot
        axs[0].plot(t, signal, color='#00CC96')
        axs[0].set_title('Original Signal')
        axs[0].set_xlabel('Time (s)')
        axs[0].set_ylabel('Amplitude')
        axs[0].grid(True)
        
        # Fourier Transform plot
        axs[1].plot(original_freqs, original_magnitude)
        axs[1].set_title('Fourier Transformation of Original Signal')
        axs[1].set_xlabel('Frequency (Hz)')
        axs[1].set_ylabel('Magnitude')
        axs[1].grid(True)
        
        plt.tight_layout()
        plt.show()

# Create the slider widget
frequency_slider = widgets.IntRangeSlider(
    value=[5, 10],
    min=1,
    max=50,
    step=1,
    description='Frequency Range (Hz)',
    style={'description_width': 'initial'},
    layout=Layout(width='90%'),
    continuous_update=False
)


def response(change):
        fig = update_plot(frequency_slider.value)

# Observe changes in the slider and update the plot accordingly
frequency_slider.observe(response, names='value')

# Style the HTML element
html_label = widgets.HTML(
    value="""
        <h2 style='font-weight: bold; font-size: 30px; text-align: center;'>Continuous Time Fourier Transform – CTFT</h2>
    """
)

# Display the slider and the outputs
vbox_layout = Layout(display='flex', flex_flow='column', align_items='center', width='100%')
ui = widgets.VBox([html_label,
              frequency_slider,
              graph_output0], layout=vbox_layout)

out = widgets.interactive_output(update_plot, {'selected_frequencies_first': frequency_slider})

clear_output(wait=True)  # Clear the previous plot
display(ui, out)
<<<<<<< HEAD <<<<<<< HEAD
=======
>>>>>>> d546939e94d695c6e9c44c085bc5ec2e2e14440e =======
>>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

By sampling \(x(t)\) at a rate \(f_s=1/T_s\), a discrete-time signal \(x(nT_s)\) is produced. Mathematically, it is represented as a series of delta functions

\[ x_\delta (t)=\sum_{n=-\infty}^{\infty}x(nT_s)\delta(t-nT_s)=x(t)\sum_{n=-\infty}^{\infty}\delta(t-nT_s) \]

with Fourier transform

\[ X_\delta (f)=\sum_{n=-\infty}^{\infty}x(nT_s)e^{-j2\pi fnT_s}=X(f)*1/T_s\sum_{n=-\infty}^{\infty}\delta(f-k/T_s)=1/T_s\sum_{n=-\infty}^{\infty}X(f-k/T_s) \]

which is a periodic function.

# Παράμετροι δειγματοληψίας
T = 1 / fs  # Sampling period
L = 1000  # Number of samples

# Create time vector
t = np.arange(0, L) * T

def plot_sampled_signal(f_range, n_samples):
    # Generate the signal with the selected frequencies
    signal = np.zeros(len(t))
    for f in range(f_range[0], f_range[1] + 1):
        signal += np.sin(2 * np.pi * f * t)
    
   # Correct downsampling approach
    sample_rate = n_samples  # Directly using selected_samples as an integer
    downsampling_factor = int(fs / sample_rate)  # Downsampling factor

    # Select every nth sample from the original time vector to match the downsampled signal
    sample_points = t[::downsampling_factor]

    # Downsample the signal by selecting every nth sample
    downsampled_signal = signal[::downsampling_factor]

    # Recalculate L for the downsampled signal if needed
    L_downsampled = len(downsampled_signal)

    # Calculate the DFT of the downsampled signal
    sampled_signal_fft = fft(downsampled_signal)

    # Calculate frequencies for the FFT of the downsampled signal
    # Note: The new sampling period is the inverse of the new sampling rate
    sampled_freqs = fftfreq(L_downsampled, 1/sample_rate)[:L_downsampled // 2]

    # Calculate magnitude of Fourier coefficients (amplitude) for the downsampled signal
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    n = len(downsampled_signal)
    T = 1 / sample_rate  # Recalculate the sampling period with the selected sample rate

    # After calculating the magnitude of Fourier coefficients
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    # Create a new frequency vector that includes the replicated frequencies
    # First, calculate the original frequency bins for the positive frequencies
    original_sampled_freqs = fftfreq(L_downsampled, T)[:L_downsampled // 2]

    shifted_freqs_negative = original_sampled_freqs - 1/T  # Shift for -1/Ts
    shifted_freqs_positive = original_sampled_freqs + 1/T  # Shift for +1/Ts

    # Concatenate the original and shifted (replicated) frequencies and magnitudes
    # This includes the original frequencies, and the replications at -1/Ts and +1/Ts
    replicated_freqs = np.concatenate([shifted_freqs_negative, original_sampled_freqs, shifted_freqs_positive])
    replicated_magnitude = np.concatenate([sampled_magnitude, sampled_magnitude, sampled_magnitude])

    # Sort the replicated frequencies and magnitudes in ascending order for proper plotting
    sorted_indices = np.argsort(replicated_freqs)
    sorted_replicated_freqs = replicated_freqs[sorted_indices]
    sorted_replicated_magnitude = replicated_magnitude[sorted_indices]
    
    # Plotting
    fig, axs = plt.subplots(1, 2, figsize=(15, 5))
    
    axs[0].scatter(sample_points, downsampled_signal, color='#00CC96', s=12)
    axs[0].set_title('Sampled Signal')
    axs[0].set_xlabel('Time (s)')
    axs[0].set_ylabel('Amplitude')
    axs[0].grid(True)
    
    axs[1].plot(sorted_replicated_freqs, sorted_replicated_magnitude, color='#1F77B4')
    axs[1].set_title('Fourier Transformation of Sampled Signal')
    axs[1].set_xlabel('Frequency (Hz)')
    axs[1].set_ylabel('Magnitude')
    axs[1].grid(True)  

    plt.tight_layout()
    plt.show()

# Widgets
frequency_slider = widgets.IntRangeSlider(
    value=[5, 10],
    min=1,
    max=50,
    step=1,
    description='Frequency Range (Hz):',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'},
    continuous_update=False
)

samples_slider = widgets.IntSlider(
    value=500,
    min=100,
    max=1000,
    step=25,
    description='Samples:',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'},
    continuous_update=False
)

html_label = widgets.HTML(
    value="""
        <h2 style='font-weight: bold; font-size: 30px; text-align: center;'>Signal Reconstruction in Digital Signal Processing</h2>
    """
)

vbox_layout = Layout(display='flex', flex_flow='column', align_items='center')
ui = widgets.VBox([html_label, frequency_slider, samples_slider], layout=vbox_layout)

out = widgets.interactive_output(plot_sampled_signal, {'f_range': frequency_slider, 'n_samples': samples_slider})

clear_output(wait=True)  # Clear the previous plot
display(ui, out)
<<<<<<< HEAD <<<<<<< HEAD
=======
>>>>>>> d546939e94d695c6e9c44c085bc5ec2e2e14440e =======
>>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

For band-limited signals \(x(t)\) with bandwidth W, assuming that the sampling rate \(fs ≥ 2W\), it holds that \(X(f) = T_s X_\delta(f)\), \(0 ≤ f ≤ W\), meaning, the signal \(X(f)\) results from passing the sampled \(x_\delta(t)\) through an ideal low-pass filter with gain \(T_s\). From the previous schema, it becomes evident that if the sampling is done with a frequency less than twice the highest frequency \(W\) of the signal (undersampling), then “images” of the spectrum from higher frequencies appear in the signal’s frequency domain, preventing the accurate restoration of the original continuous-time signal. This phenomenon is called aliasing, and the error during the restoration of the original signal is referred to as aliasing error. Sampling in the time domain is the basis for defining the Discrete Time Fourier Transform (DTFT). For a series of discrete numbers \(x[n]\), the discrete-time Fourier transform is defined as:

\[ X_d(\phi) \triangleq \sum_{n=-\infty}^{\infty} x[n] \exp (-j 2 \pi n \phi) \]

The DTFT is a periodic function with a period of \(1\), therefore, it is sufficient to compute it in the frequency interval \([0,1]\) or equivalently \([-½,½]\). Note that the DTFT, although it arises from a series of discrete numbers \(x[n]\), is a continuous function of the variable \(\phi\) as illustrated in the next schema.

# Παράμετροι δειγματοληψίας
T = 1 / fs  # Sampling period
L = 1000  # Number of samples

# Create time vector
t = np.arange(0, L) * T

def plot_sampled_signal(selected_frequencies, selected_samples):
    # Generate the signal with the selected frequencies
    signal = np.zeros(len(t))
    for f in range(selected_frequencies[0], selected_frequencies[1] + 1):
        signal += np.sin(2 * np.pi * f * t)
    
   # Correct downsampling approach
    sample_rate = selected_samples  # Directly using selected_samples as an integer
    downsampling_factor = int(fs / sample_rate)  # Downsampling factor

    # Downsample the signal by selecting every nth sample
    downsampled_signal = signal[::downsampling_factor]

    # Recalculate L for the downsampled signal if needed
    L_downsampled = len(downsampled_signal)

    # Calculate the DFT of the downsampled signal
    sampled_signal_fft = fft(downsampled_signal)

    # Calculate magnitude of Fourier coefficients (amplitude) for the downsampled signal
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    n = len(downsampled_signal)
    T = 1 / sample_rate  # Recalculate the sampling period with the selected sample rate

    # After calculating the magnitude of Fourier coefficients
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    # Create a new frequency vector that includes the replicated frequencies
    # First, calculate the original frequency bins for the positive frequencies
    original_sampled_freqs = fftfreq(L_downsampled, T)[:L_downsampled // 2]

    shifted_freqs_negative = original_sampled_freqs - 1/T  # Shift for -1/Ts
    shifted_freqs_positive = original_sampled_freqs + 1/T  # Shift for +1/Ts

    # Concatenate the original and shifted (replicated) frequencies and magnitudes
    # This includes the original frequencies, and the replications at -1/Ts and +1/Ts
    replicated_freqs = np.concatenate([shifted_freqs_negative, original_sampled_freqs, shifted_freqs_positive])
    replicated_magnitude = np.concatenate([sampled_magnitude, sampled_magnitude, sampled_magnitude])

    # Sort the replicated frequencies and magnitudes in ascending order for proper plotting
    sorted_indices = np.argsort(replicated_freqs)
    sorted_replicated_freqs = replicated_freqs[sorted_indices]
    sorted_replicated_magnitude = replicated_magnitude[sorted_indices]

    # Generate an array of sample indices for the downsampled signal
    sample_indices = np.arange(len(downsampled_signal))

    # Normalize the frequency values
    normalized_freqs = original_sampled_freqs / (1/T)

    # Since you previously concatenated and sorted for replication, ensure to apply normalization there as well
    normalized_replicated_freqs = np.concatenate([
        (shifted_freqs_negative / (1/T)),  # Normalize -1/Ts shifted frequencies
        normalized_freqs,  # Already normalized original frequencies
        (shifted_freqs_positive / (1/T))   # Normalize +1/Ts shifted frequencies
    ])

    # Sort the normalized and replicated frequencies for proper plotting
    sorted_indices = np.argsort(normalized_replicated_freqs)
    sorted_normalized_replicated_freqs = normalized_replicated_freqs[sorted_indices]
    sorted_replicated_magnitude = replicated_magnitude[sorted_indices]

    # Create the plots
    fig, axs = plt.subplots(1, 2, figsize=(15, 5))
    
    # Plot the downsampled signal
    axs[0].scatter(sample_indices, downsampled_signal, color='#00CC96', s=12)
    axs[0].set_title('Sampled Signal')
    axs[0].set_xlabel('Sample Number')
    axs[0].set_ylabel('Amplitude')
    axs[0].grid(True)

    # Plot the Fourier Transform of the downsampled signal
    axs[1].plot(sorted_normalized_replicated_freqs, sorted_replicated_magnitude, color='#1F77B4')
    axs[1].set_title('Normalized Fourier Transformation of Sampled Signal')
    axs[1].set_xlabel('Normalized Frequency (f/fs)')
    axs[1].set_ylabel('Magnitude')
    axs[1].grid(True)
    
    plt.tight_layout()
    plt.show()

# IPyWidgets sliders
frequency_slider = widgets.IntRangeSlider(
    value=[5, 10],
    min=1,
    max=50,
    step=1,
    description='Frequency range (Hz):',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'},
    continuous_update=False
)

samples_slider = widgets.IntSlider(
    value=500,
    min=100,
    max=1000,
    step=25,
    description='Samples:',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'},
    continuous_update=False
)

html_label = widgets.HTML(
    value="""
        <h2 style='font-weight: bold; font-size: 30px; text-align: center;'>Discrete Time Fourier Transform – DTFT</h2>
    """
)

# Display the sliders and output
vbox_layout = Layout(display='flex', flex_flow='column', align_items='center')
ui = widgets.VBox([html_label, frequency_slider, samples_slider], layout=vbox_layout)

# Interaction between widgets and function
out = widgets.interactive_output(
    plot_sampled_signal, 
    {'selected_frequencies': frequency_slider, 'selected_samples': samples_slider}
)

clear_output(wait=True)  # Clear the previous plot
display(ui, out)
<<<<<<< HEAD <<<<<<< HEAD
=======
>>>>>>> d546939e94d695c6e9c44c085bc5ec2e2e14440e ======= >>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

With the series of discrete numbers arising as a result of sampling, \(x[n]=x(nT_s)\), the DTFT and the Fourier transform \(X_\delta(f)\) of the sampled signal are connected through the correspondence \(\phi ↔ f/f_s\). The common practice is to represent the ratio \(f/f_s\) as normalized frequency \(\phi\) (\(f_D\), in your notes) and the real frequencies to arise as multiples of it (usually fractional). To link the DTFT with the Fourier transform \(X(f)\) of the signal, an additional conversion to the sampling period by multiplying by \(T_s\) (or dividing by \(f_s\)) is required. Analogous to sampling signals in time, we can sample in the frequency domain by taking discrete values \(X(kf_o)\) of the Fourier transform corresponding to a frequency analysis \(f_o=1/T_o\). This is equivalent to the periodic repetition of the continuous-time signal \(x(t)\) every \(Τ_ο\), since the periodic signal

\[ x_p(t)=\sum_{n=-\infty}^{\infty} x\left(t-n T_{o}\right) \]

has Fourier transform

\[ X(f) \sum_{n=-\infty}^{\infty} \exp \left(-j 2 \pi f n T_{o}\right)=X(f) \frac{1}{T_{o}} \sum_{k=-\infty}^{\infty} \delta\left(f-k / T_{o}\right)=\frac{1}{T_{o}} \sum_{k=-\infty}^{\infty} X\left(k / T_{o}\right) \delta\left(f-k / T_{o}\right) \]

Therefore, \(X[k] = X(kf_o)/Τ_o\) are the coefficients of the Fourier series expansion of the periodic signal \(x_p(t)\). Obviously, for signals \(x(t)\) of finite duration, where \(x(t)=0\) for \(|t| ≥ T\), assuming that the period \(T_o ≥ 2T\), it holds that \(x(t) = x_p(t)\) for \(|t| ≤ T\). In practice, signals have a very long duration to be analyzed in their entirety. Thus, we apply a rectangular time window to retain only their most significant part for the observation period and \(x(t)= 0\) elsewhere. In computing the DTFT \(X_d(\phi)\) of such a truncated signal, instead of an infinite sum, we are limited to a finite-length \(L\) series of numbers \(x[n]\), therefore

\[ X_d(\phi)=\sum_{n=0}^{L-1} x[n] \exp (-j 2 \pi n \phi) \]

The sampling of \(X_d(\phi)\) in the frequency domain at \(N\) equidistant normalized frequencies \(0\), \(1/N\), \(2/N\), \(…\), \((N-1)/N\), gives

\[ X[k]=X_{d}\left(\frac{k}{N}\right)=\sum_{n=0}^{N-1} x[n] \exp \left(-j 2 \pi n \frac{k}{N}\right), \quad 0 \leq k \leq N-1 \]

where, if \(N≥L\), we set \(x[n]=0\) for \(n≥L\). The last relationship is recognized as the Discrete Fourier Transform (DFT), which for a finite series \(xn\), \(n=0\), \(1\), \(…\), \(N-1\), is defined as:

\[ X_{k} \triangleq \sum_{n=0}^{N-1} x_{n} \exp \left(-j 2 \pi n \frac{k}{N}\right), \quad 0 \leq k \leq N-1 \]

and its inverse is

\[ x_{n}=\frac{1}{N} \sum_{k=0}^{N-1} X_{k} \exp \left(j 2 \pi n \frac{k}{N}\right), \quad 0 \leq n \leq N-1 \]

The \(X_d(\phi)\) as DTFT is a periodic function and if the original series xn was periodic (and we didn’t apply the window), then \(X_d(\phi)\) would be zero everywhere except at the sampling points \(k/N\). That is, if we consider a finite-length series of numbers that repeats periodically, the discrete-time Fourier transform of it (DTFT) is also periodic and discrete. Moreover, the DFT and its inverse IDFT, if we did not limit the indices \(n\) and \(k\) between \(0\) and \(N-1\), would be periodic functions. Therefore, the finite series xn can be considered as a periodic discrete-time signal observed only during a period and the DFT, the series \(X_k\), as the samples with resolution \(1/N\) of the DTFT \(X_d(\phi)\) in the domain of normalized frequencies \([0,1]\), as shown in the next schema.

# Παράμετροι δειγματοληψίας
T = 1 / fs  # Sampling period
L = 1000  # Number of samples

# Create time vector
t = np.arange(0, L) * T

def plot_signals(f_range, samples, selected_N):
    # Generate the signal with the selected frequencies
    signal = np.zeros(len(t))
    for f in range(f_range[0], f_range[1] + 1):
        signal += np.sin(2 * np.pi * f * t)
    
   # Correct downsampling approach
    sample_rate = samples  # Directly using samples as an integer
    downsampling_factor = int(fs / sample_rate)  # Downsampling factor

    # Downsample the signal by selecting every nth sample
    downsampled_signal = signal[::downsampling_factor]

    # Recalculate L for the downsampled signal if needed
    L_downsampled = len(downsampled_signal)

    # Calculate the DFT of the downsampled signal
    sampled_signal_fft = fft(downsampled_signal)

    # Calculate magnitude of Fourier coefficients (amplitude) for the downsampled signal
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    n = len(downsampled_signal)
    T = 1 / sample_rate  # Recalculate the sampling period with the selected sample rate

    # After calculating the magnitude of Fourier coefficients
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    # Create a new frequency vector that includes the replicated frequencies
    # First, calculate the original frequency bins for the positive frequencies
    original_sampled_freqs = fftfreq(L_downsampled, T)[:L_downsampled // 2]

    shifted_freqs_negative = original_sampled_freqs - 1/T  # Shift for -1/Ts
    shifted_freqs_positive = original_sampled_freqs + 1/T  # Shift for +1/Ts

    # Concatenate the original and shifted (replicated) frequencies and magnitudes
    # This includes the original frequencies, and the replications at -1/Ts and +1/Ts
    replicated_freqs = np.concatenate([shifted_freqs_negative, original_sampled_freqs, shifted_freqs_positive])
    replicated_magnitude = np.concatenate([sampled_magnitude, sampled_magnitude, sampled_magnitude])

    # Sort the replicated frequencies and magnitudes in ascending order for proper plotting
    sorted_indices = np.argsort(replicated_freqs)
    sorted_replicated_freqs = replicated_freqs[sorted_indices]
    sorted_replicated_magnitude = replicated_magnitude[sorted_indices]

    # Generate an array of sample indices for the downsampled signal
    sample_indices = np.arange(len(downsampled_signal))

    # Normalize the frequency values
    normalized_freqs = original_sampled_freqs / (1/T)

    # Since you previously concatenated and sorted for replication, ensure to apply normalization there as well
    normalized_replicated_freqs = np.concatenate([
        (shifted_freqs_negative / (1/T)),  # Normalize -1/Ts shifted frequencies
        normalized_freqs,  # Already normalized original frequencies
        (shifted_freqs_positive / (1/T))   # Normalize +1/Ts shifted frequencies
    ])

    # Sort the normalized and replicated frequencies for proper plotting
    sorted_indices = np.argsort(normalized_replicated_freqs)
    sorted_normalized_replicated_freqs = normalized_replicated_freqs[sorted_indices]
    sorted_replicated_magnitude = replicated_magnitude[sorted_indices]

    # Adjust the sampling interval for selecting frequencies and magnitudes
    sampled_indices = np.arange(0, len(sorted_normalized_replicated_freqs), selected_N)
    sampled_normalized_replicated_freqs = sorted_normalized_replicated_freqs[sampled_indices]
    sampled_replicated_magnitude = sorted_replicated_magnitude[sampled_indices]
    
    # Create the plots
    fig, axs = plt.subplots(1, 2, figsize=(15, 5))

    # Time-domain signal plot
    axs[0].scatter(sample_indices, downsampled_signal, color='#00CC96', s=12)
    axs[0].set_title('Sampled Signal')
    axs[0].set_xlabel('Sample Number')
    axs[0].set_ylabel('Amplitude')
    axs[0].grid(True)

    
    axs[1].stem(sampled_normalized_replicated_freqs, sampled_replicated_magnitude, linefmt='#1F77B4', markerfmt='o', basefmt=" ")
    axs[1].set_title('Sampled Normalized Fourier Transformation')
    axs[1].set_xlabel('Normalized Frequency (f/fs)')
    axs[1].set_ylabel('Magnitude')
    axs[1].grid(True)

    plt.tight_layout()
    plt.show()

# Define interactive widgets
frequency_slider = widgets.IntRangeSlider(
    value=[5, 40],
    min=1,
    max=51,
    step=1,
    description='Frequency range (Hz):',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'}
)

samples_slider = widgets.IntSlider(
    value=500,
    min=100,
    max=1000,
    step=25,
    description='Samples:',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'}
)

N_slider = widgets.IntSlider(
    value=1,
    min=1,
    max=25,
    step=1,
    description='N:',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'}
)

html_label = widgets.HTML(
    value="""
        <h2 style='font-weight: bold; font-size: 30px; text-align: center;'>Properties of the Discrete Fourier Transform in Digital Signal Processing</h2>
    """
)

# Display the sliders and output
vbox_layout = Layout(display='flex', flex_flow='column', align_items='center')
ui = widgets.VBox([html_label, frequency_slider, samples_slider, N_slider], layout=vbox_layout)
out = widgets.interactive_output(plot_signals, {'f_range': frequency_slider, 'samples': samples_slider, 'selected_N': N_slider})

# Display the widgets and output
clear_output(wait=True)  # Clear the previous plot
display(ui, out)
<<<<<<< HEAD <<<<<<< HEAD
=======
>>>>>>> d546939e94d695c6e9c44c085bc5ec2e2e14440e ======= >>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

Spectral Analysis#

To calculate the energy or power of the waveform \(x(t)\), depending on the case of the signal, it holds

\[ E_{X}=\int_{-\infty}^{\infty} x^{2}(t) d t=\int_{-\infty}^{\infty}|X(f)|^{2} d f \]
\[ P_{X}=\lim _{T \rightarrow \infty} \frac{1}{T} \int_{-T / 2}^{T / 2} x^{2}(t) d t=\int_{-\infty}^{\infty} S_{X}(f) d f \]

where for power signals \(S_Χ(f)\) is the power spectral density (PSD) of \(x(t)\). For discrete-time signals resulting from sampling of \(x(t)\) with period \(T_s\), the corresponding relations for calculating the energy or power become

\[ E_{X}=T_{s} \sum_{n=-\infty}^{\infty} x^{2}[n] \]
\[ P_{X}=\lim _{N \rightarrow \infty} \frac{1}{2 N+1} \sum_{n=-N}^{N} x^{2}[n] \]

A simple way to estimate the power spectral density of the waveform \(x(t)\) is to take the DTFT of the signal’s samples and then square the magnitude of the result. This estimator is called a periodogram. The periodogram of a finite-length \(L\) signal \(x[n]\) is defined as

\[ P_{x x}(f) \triangleq \frac{\left|X_{d}\left(f / f_{s}\right)\right|^{2}}{f_{s} L} \]

where \(X_d(\phi)\) the DTFT of the signal. As the length \(L\) tends to infinity, the periodogram \(P_{xx}(f)\) approaches the power spectral density \(S_Χ(f)\). The calculation of the periodogram at a finite number of frequencies \(kf_s/N\), \(k=0\), \(1\), \(…\), \(N\) gives

\[ P_{x x}[k]=\frac{\left|X_{k}\right|^{2}}{f_{s} L}, \quad k=0,1, \ldots, N-1 \]

where \(X_k\) is the DFT of the finite-length \(L\) series of signal samples. The power of the signal is then

\[ P_{X}=\frac{1}{f_{s} L} \sum_{k=0}^{N-1}\left|X_{k}\right|^{2} f_{o}=\frac{1}{N L} \sum_{k=0}^{N-1}\left|X_{k}\right|^{2}=\frac{1}{L} \sum_{n=0}^{L-1}\left|x_{n}\right|^{2} \]

where the last equality results from the Parseval’s theorem, which for the DFT case is expressed as:

\[ \sum_{n=0}^{N-1}\left|x_{n}\right|^{2}=\frac{1}{N} \sum_{k=0}^{N-1}\left|X_{k}\right|^{2} \]

In the special case of periodic signals we have

\[ S_x(f) = \sum_{k=-\infty}^{\infty} |X[k]|^2 \delta(f - k/T_0) \]
\[ P_x = \sum_{k=-\infty}^{\infty} |X[k]|^2 \]

where \(X[k]\) are the coefficients of the Fourier series expansion and \(T_o\) the period of the signal.

Application in Practice#

Create a sine wave signal#

Code
import numpy as np
import matplotlib.pyplot as plt

# ==============================================================================
# 2.1 Create a Sinusoidal Signal
# ==============================================================================
Fs = 2000                  # Sampling frequency in Hz
Ts = 1 / Fs                # Sampling period in seconds
T = 0.1                    # Signal duration in seconds
t = np.arange(0, T, Ts)    # Time vector for signal
A = 1                      # Signal amplitude
x = A * np.sin(2 * np.pi * 100 * t)  # Generate sinusoidal signal
L = len(x)                 # Length of the signal

# Plot the sinusoidal signal in time domain
plt.figure()
plt.plot(t, x)
plt.title('Sinusoidal Signal')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.show(block=False)
plt.pause(1)  # Pause to view the plot, press any key in the console to continue




% ==============================================================================
% 2.1 Create a Sinusoidal Signal
% ==============================================================================
Fs = 2000;                 % Sampling frequency in Hz
Ts = 1 / Fs;               % Sampling period in seconds
T = 0.1;                   % Signal duration in seconds
t = 0:Ts:T-Ts;             % Time vector for signal
A = 1;                     % Signal amplitude
x = A * sin(2 * pi * 100 * t); % Generate sinusoidal signal
L = length(x);             % Length of the signal
plot(t, x);                % Plot the sinusoidal signal in time domain
title('Sinusoidal Signal');% Title for the plot
xlabel('Time (s)');        % X-axis label
ylabel('Amplitude');       % Y-axis label
pause;                     % Pause to view the plot, press any key to continue
../_images/162abd163ab53285efdeebfd0b8c1d8107ceacc32264c47e38badc21475160ff.png

Plot the discrete Fourier transform of the sine wave signal#

Code
# ==============================================================================
# 2.2 Plot Fourier Transform (FT) of the Signal
# ==============================================================================
N = 1 * L                  # Length of Fourier Transform
Fo = Fs / N                # Frequency resolution
Fx = np.fft.fft(x, N)      # Discrete Fourier Transform (DFT) of the signal
freq = np.arange(0, N) * Fo  # Frequency vector

# Plot the magnitude of the DFT
plt.figure()
plt.plot(freq, np.abs(Fx))
plt.title('FFT of the Signal')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.show(block=False)
plt.pause(1)  # Pause to view the plot, press any key in the console to continue
plt.axis([0, 100, 0, L/2])
plt.pause(1)



% ==============================================================================
% 2.2 Plot Fourier Transform (FT) of the Signal
% ==============================================================================
N = 1 * L;                 % Length of Fourier Transform
Fo = Fs / N;               % Frequency resolution
Fx = fft(x, N);            % Discrete Fourier Transform (DFT) of the signal
freq = (0:N-1) * Fo;       % Frequency vector
plot(freq, abs(Fx));       % Plot the magnitude of the DFT
title('FFT of the Signal');% Title for the plot
xlabel('Frequency (Hz)');  % X-axis label
ylabel('Magnitude');       % Y-axis label
pause;                     % Pause to view the plot, press any key to continue
axis([0 100 0 L/2]);       % Set axis limits
pause;                     % Pause, press any key to continue
../_images/eddb197c3f88cd80ec8c06b0a19f4875d513359bbb1c07dc6e0663b07f5ff10a.png

Plot the periodogram#

Code
# ==============================================================================
# 2.3 Plot Signal Periodogram
# ==============================================================================
power = (Fx * np.conj(Fx)) / (Fs * L)  # Calculate spectral density

plt.figure()
plt.plot(freq, power)
plt.title('Periodogram')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.show(block=False)
plt.pause(1)




% ==============================================================================
% 2.3 Plot Signal Periodogram
% ==============================================================================
power = Fx .* conj(Fx) / Fs / L; % Calculate spectral density
plot(freq, power);         % Plot the periodogram
title('{\bf Periodogram}');% Title for the plot with bold characters
xlabel('Frequency (Hz)');  % X-axis label
ylabel('Power');           % Y-axis label
../_images/53bfce46baf9225a64949aee8e46dcfd86b3e2b0252db55eacc5ee4c7d8fbe35.png

Calculate the power of the sine wave signal#

Code
# ==============================================================================
# 2.4 Calculate Signal Power
# ==============================================================================
power_theory = A**2 / 2                # Theoretical power based on signal amplitude
dB = 10 * np.log10(power_theory)       # Convert power to decibels (dB)
power_time_domain = np.sum(np.abs(x)**2) / L  # Calculate power in time domain
power_frequency_domain = np.sum(power) * Fo  # Calculate power in frequency domain

# Display calculated power values
print(f'Power (Theory): {power_theory}')
print(f'Power (dB): {dB}')
print(f'Power (Time Domain): {power_time_domain}')
print(f'Power (Frequency Domain): {power_frequency_domain}')


% ==============================================================================
% 2.4 Calculate Signal Power
% ==============================================================================
power_theory = A^2 / 2;            % Theoretical power based on signal amplitude
dB = 10 * log10(power_theory);     % Convert power to decibels (dB)
power_time_domain = sum(abs(x).^2) / L; % Calculate power in time domain
power_frequency_domain = sum(power) * Fo; % Calculate power in frequency domain

% Display calculated power values
disp(['Power (Theory): ', num2str(power_theory)]);
disp(['Power (dB): ', num2str(dB)]);
disp(['Power (Time Domain): ', num2str(power_time_domain)]);
disp(['Power (Frequency Domain): ', num2str(power_frequency_domain)]);

Part 3: Application A#

Step 1: Signal creation and display#

Next, you will apply what you have learned in a more complex example of signal generation that includes modulation and noise addition. For convenience, incomplete PYTHON/MATLAB code with comments is provided, which you need to complete and save.

Code
    # ==============================================================================
    # Part 1: Create the signal
    # ==============================================================================
    import numpy as np
    import matplotlib.pyplot as plt

    # Close all figure windows
    ______

    # Clear workspace
    # In Python, variables need to be explicitly deleted if needed. Typically, just overwrite or ignore.

    # Clear command window
    # This is not applicable in Python in the same way as in MATLAB.

    Fs = 4000  # Sampling frequency in Hz
    Ts = ___  # Sampling period
    L = 3000  # Signal length (number of samples)
    T = L * Ts  # Signal duration
    t = np.arange(0, L) * Ts  # Timestamps of signal calculation

    # Signal creation
    x = np.sin(2 * np.pi * 200 * t) \
        + 0.3 * np.sin(2 * np.pi * 300 * (t - 2)) \
        + np.sin(2 * np.pi * 400 * t)

    # Plot signal x in time
    plt.figure(1)
    plt.plot(t, x)
    plt.title('Time domain plot of x')
    plt.xlabel('t (sec)')
    plt.ylabel('Amplitude')
    plt.show(block=False)  # block=False allows the script to continue while the plot is shown
    plt.pause(1)  # Pause for 1 second; adjust as needed

    plt.axis([0, 0.3, -2, 2])  # Adjust axes
    plt.pause(1)  # Pause for 1 second; adjust as needed

    # Calculate Fourier Transform
    N = 2 ** np.ceil(np.log2(L)).astype(int)  # FFT length
    Fo = ___  # Frequency analysis step
    f = np.arange(0, N) * Fo  # Frequency vector
    X = ___  # DFT for N points, fill in with the correct function

    # Plot signal in frequency domain
    plt.figure(2)
    plt.plot(f[:N // 2], np.abs(X[:N // 2]))  # Plot signal only for positive frequency values
    plt.title('Frequency domain plot of x')
    plt.xlabel('f (Hz)')
    plt.ylabel('Amplitude')
    plt.show(block=False)
    plt.pause(1)

    # For a two-sided signal plot
    plt.figure(3)
    f = f - Fs / 2  # Shift frequency vector
    X = np.fft.fftshift(X)  # Shift the zero to the center
    plt.plot(f, np.abs(X))
    plt.title('Two sided spectrum of x')
    plt.xlabel('f (Hz)')
    plt.ylabel('Amplitude')
    plt.show(block=False)
    plt.pause(1)

    # Calculate power of signal
    power = (X * np.conj(X)) / N / L
    plt.figure(4)
    plt.plot(f, power)
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Power')
    plt.title('Periodogram')
    plt.show(block=False)
    plt.pause(1)

% ==============================================================================
% Part 1 Create the signal
% ==============================================================================
_________ % close all figure windows
_________ % clear workspace
_________ % clear command window
Fs=4000; % sampling frequency 4000 Hz
Ts=____; % sampling period
L=3000; % signal length (number of samples)
T=L*Ts; % signal duration
t=0:Ts:(L-1)*Ts; % timestamps of signal calculation
x=sin(2*pi*200*t)... % sinusoidal signal of frequency 200 Hz
+ 0.3*sin(2*pi*300*(t-2))... % plus sinusoidal signal of frequency 300 Hz
+ sin(2*pi*400*t); % plus sinusoidal signal of frequency 400 Hz
% Plot signal x in time
figure(1) % open a new figure window
plot(t,x) % plot the signal
title('Time domain plot of x') % plot title
xlabel('t (sec)') % x axis label
ylabel('Amplitude') % y axis label
pause % pause, press any key to continue
axis([0 0.3 -2 2]) % x axis values from 0 to 0.3 sec and
% y axis values from -2 to 2
pause % pause, press any key to continue
% Calculate Fourier Transform
N = 2^nextpow2(L); % FFT length
% nextpow2 calculates next power of 2
% grater or equal to L
Fo=____; % frequency analysis
f=(0:N-1)*Fo; % frequency vector
X=________; % DFT for Ν points
% Plot signal in frequency domain
% As signal is real you can plot the one sided signal(positive frequency values)
figure(2) % open a new figure window
plot(f(1:_____)),abs(X(1:_____))) % plot signal only for positive
%frequency values
title('Frequency domain plot of x') % plot title
xlabel('f (Hz)') % set x axis label
ylabel('Amplitude') % set y axis label
pause % pause, press any key to continue
% in order to plot the 2-sided signal you can use fftshift matlab
% function, type help fftshift for more details
figure(3) % open a new figure window
f=f-Fs/2; % shift frequency vector to the left by –Fs/2
X=fftshift(X); % create 2-sided signal
plot(f,abs(X));title('Two sided spectrum of x'); xlabel('f (Hz)');
ylabel('Amplitude')
pause % pause, press any key to continue
% Calculate power of signal
power=X.*conj(X)/N/L; % Calculate power of signal
figure(4) % open a new figure window
plot(f,power) % plot signal power to frequency
xlabel('Frequency (Hz)') % x axis label
ylabel('Power') % y axis label
title('{\bf Periodogram}') % plot title
pause

# Callback function to update graphs
def update_graph(selected_frequencies_part1, selected_Fo):
    # Part 1: Create the signal
    Fs = selected_frequencies_part1        # Sampling frequency 1000 Hz
    Ts = 1 / Fs                      # Sampling period
    L = 1000                         # Length of signal (number of samples)
    T = L * Ts                       # Duration of signal
    t = np.arange(0, (L - 1) * Ts, Ts)  # Time vector

    global new_x 
    new_x = np.sin(2 * np.pi * (selected_Fo-30) * t) \
        + 0.8 * np.sin(2 * np.pi * (selected_Fo+20) * (t - 2)) \
        + np.sin(2 * np.pi * selected_Fo * t)  # 60 Hz component
    
    # Plotting
    fig, axs = plt.subplots(4, 1, figsize=(12, 20))
    
    # Time domain plot
    axs[0].plot(t, new_x, color='#00CC96')
    axs[0].set_title('Time domain plot of x')
    axs[0].set_xlabel('t (sec)')
    axs[0].set_ylabel('Amplitude')
    axs[0].grid(True)
    
    # Fourier transform
    def nextpow2(i):
        n = 1
        while n < i:
            n *= 2
        return n

    N = nextpow2(L)                 # Length of Fourier transform
    Fo = Fs / N                     # Frequency resolution
    f = np.arange(0, N) * Fo        # Frequency vector
    X = np.fft.fft(new_x, N)            # Compute DFT for N points

    # Frequency domain plot
    axs[1].plot(f[1:N], abs(X[1:N]), color='#1F77B4')
    axs[1].set_title('Frequency domain plot of x')
    axs[1].set_xlabel('f (Hz)')
    axs[1].set_ylabel('Amplitude')
    axs[1].grid(True)

    # Shift frequencies to center
    f = f - Fs / 2
    X = np.fft.fftshift(X)

    # Two-sided spectrum of x
    f_shifted = f 

    # Two-sided spectrum plot
    axs[2].plot(f_shifted, abs(X), color='#1F77B4')
    axs[2].set_title('Two sided spectrum of x')
    axs[2].set_xlabel('f (Hz)')
    axs[2].set_ylabel('Amplitude')
    axs[2].grid(True)

    # Calculate power
    power = np.multiply(X, np.conj(X)) / N / L

    # Periodogram plot
    axs[3].plot(f_shifted, power.real, color='#1F77B4')
    axs[3].set_title('Periodogram')
    axs[3].set_xlabel('Frequency (Hz)')
    axs[3].set_ylabel('Power')
    axs[3].grid(True)

    plt.tight_layout()
    plt.show()

# Create interactive widgets
single_frequency_slider = widgets.IntSlider(
    min=100,
    max=2000,
    step=100,
    value=1000,
    description='Sampling Frequency (Fs):',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'},
    continuous_update=False
)

Fo_slider = widgets.IntSlider(
    min=40,
    max=400,
    step=10,
    value=100,
    description='Frequency (Fo):',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'},
    continuous_update=False
)

html_label = widgets.HTML(
    value="""

    """
)

# Display the sliders and output
vbox_layout = Layout(display='flex', flex_flow='column', align_items='center')

ui = widgets.VBox([html_label, single_frequency_slider, Fo_slider], layout=vbox_layout)
out = widgets.interactive_output(update_graph, {'selected_frequencies_part1': single_frequency_slider, 'selected_Fo': Fo_slider})

# Display the widgets and output
clear_output(wait=True)  # Clear the previous plot
display(ui, out)
<<<<<<< HEAD <<<<<<< HEAD
=======
>>>>>>> d546939e94d695c6e9c44c085bc5ec2e2e14440e ======= >>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

Step 2: Add noise to the signal#

Complete the code to create the noise signal n using the randn function. The noise vector n should be the same size as the sine wave x from the first part. Plot the noise signal in the interval from 0 to 0.2 sec and scale from -2 to 2. Calculate the periodogram of n and plot the power spectral density of the noise signal. Add the noise signal and x to get the noisy signal s. Plot the noisy signal s in the time domain in the area from 0 to 0.2 sec and scale from -2 to 2 as well as its bilateral spectrum.

Fs=1000                    # συχνότητα δειγματοληψίας 1000 Hz
Ts=1/Fs                    # περίοδος δειγματοληψίας
L=1000                     # μήκος σήματος (αριθμός δειγμάτων)
T=L*Ts                     # διάρκεια σήματος
t=np.arange(0,(L-1)*Ts,Ts) # χρονικές στιγμές υπολογισμού του σήματος

def nextpow2(i):
    # Compute the next highest power of 2
    n = 1
    while n < i: n *= 2
    return n
             

# Function to update plots
def update_plots(selected_frequencies_part2):
    new_x=np.sin(2*np.pi*(selected_frequencies_part2-30)*t) + 0.8*np.sin(2*np.pi*(selected_frequencies_part2+20)*(t-2))+ np.sin(2*np.pi*(selected_frequencies_part2)*t);         
    rand_n = np.random.randn(np.size(new_x))
    
    # Plotting
    fig, axs = plt.subplots(4, 1, figsize=(12, 20))

    # Time domain plot of n
    axs[0].plot(t, rand_n, color='#00CC96')
    axs[0].set_title('Time domain plot of n')
    axs[0].set_xlabel('t (sec)')
    axs[0].set_ylabel('Amplitude')
    axs[0].grid(True)
    
    # Correction for N calculation using bitwise operator
    N = 2^nextpow2(L)
    Fo = Fs / N   
    f = (np.arange(0, N)) * Fo
    f_shifted = f - Fs/2
    rand_N = np.fft.fft(rand_n, N)
    rand_N = np.fft.fftshift(rand_N)
    power_n = np.multiply(rand_N, np.conj(rand_N)) / N / L
    # Frequency domain plot of x
    axs[1].plot(f_shifted, power_n.real, color='#1F77B4')
    axs[1].set_title('Frequency domain plot of x')
    axs[1].set_xlabel('f (Hz)')
    axs[1].set_ylabel('Amplitude')
    axs[1].grid(True)
    
    # Two sided spectrum of x
    s = new_x + rand_n
    axs[2].plot(t, s, color='#00CC96')
    axs[2].set_title('Two sided spectrum of x')
    axs[2].set_xlabel('t (sec)')
    axs[2].set_ylabel('Amplitude')
    axs[2].grid(True)
    
    # Two sided spectrum of s
    S = np.fft.fft(s, N)
    S = np.fft.fftshift(S)
    axs[3].plot(f_shifted, np.abs(S), color='#1F77B4')
    axs[3].set_title('Two sided spectrum of s')
    axs[3].set_xlabel('f (Hz)')
    axs[3].set_ylabel('Magnitude')
    axs[3].grid(True)
    
    plt.tight_layout()
    plt.show()

# Create the slider widget
fs_slider = widgets.IntSlider(
    value=1000,
    min=100,
    max=2000,
    step=100,
    description='Sampling Frequency (Hz):',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'},
    continuous_update=False
)

html_label = widgets.HTML(
    value="""

    """
)

# Display the sliders and output
vbox_layout = Layout(display='flex', flex_flow='column', align_items='center')

ui = widgets.VBox([html_label, fs_slider], layout=vbox_layout)
out = widgets.interactive_output(update_plots, {'selected_frequencies_part2': fs_slider})

# Display the slider
clear_output(wait=True)  # Clear the previous plot
display(ui,out)
<<<<<<< HEAD <<<<<<< HEAD
=======
>>>>>>> d546939e94d695c6e9c44c085bc5ec2e2e14440e ======= >>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

Step 3: Multiplication of signals#

Complete the code for creating a sine wave signal of 100 Hz frequency and multiply it with the previous signal s. The two signals should be of the same size. Plot the result in the time domain in the area from 0 to 0.2 sec and scale from -2 to 2 as well as in the frequency domain using the fftshift function.

# Part 3. Πολλαπλασιασμός σημάτων

# Συμπληρώστε τον κώδικα δημιουργίας ενός ημιτονοειδούς σήματος συχνότητας
# 100 Hz και πολλαπλασιάστε με το προηγούμενο σήμα s.
# Τα δύο σήματα θα πρέπει να είναι του ίδιου μεγέθους.
# Σχεδιάστε το αποτέλεσμα στο πεδίο του χρόνου στην περιοχή 0 έως 0.2 sec
# και κλίμακα από -2 έως 2 καθώς και στο πεδίο της συχνότητας
# χρησιμοποιώντας τη συνάρτηση fftshift.

def update_plots(selected_frequencies_part3):
    Fο=selected_frequencies_part3
    z=np.sin(2*np.pi*Fο*t)
    new_x = np.sin(2 * np.pi * 30 * t) \
        + 0.8 * np.sin(2 * np.pi * 80 * (t - 2)) \
        + np.sin(2 * np.pi * 60 * t)  # 60 Hz component
    L = 1000  # Length of signal
    rand_n = np.random.normal(0, 1, np.size(new_x))  # Example random noise
    s = new_x + rand_n
    y= np.multiply(z,s)

    # Plotting
    fig, axs = plt.subplots(2, 1, figsize=(12, 20))

    # Time domain plot
    axs[0].plot(t, y, color='#00CC96')
    axs[0].set_title('Time domain plot of y')
    axs[0].set_xlabel('t (sec)')
    axs[0].set_ylabel('Amplitude')
    axs[0].grid(True)
    
    # Fourier transform
    def nextpow2(i):
        n = 1
        while n < i:
            n *= 2
        return n
            
    N = 2^nextpow2(L)
    Fo = Fs / N   
    f = (np.arange(0, N)) * Fo
    Y = np.fft.fft(y, N)
    f=f-Fs/2   
    Y = np.fft.fftshift(Y)

    axs[1].plot(f, np.abs(Y), color='#1F77B4')
    axs[1].set_title('Frequency domain plot of y')
    axs[1].set_xlabel('f (Hz)')
    axs[1].set_ylabel('Amplitude')
    axs[1].grid(True)

    plt.tight_layout()
    plt.show()

# Create the slider widget for Fo
Fo_slider_step3 = widgets.IntSlider(
    value=100,
    min=10,
    max=200,
    step=10,
    description='Fo (Hz):',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'},
    continuous_update=False
)

html_label = widgets.HTML(
    value="""

    """
)

# Display the sliders and output
vbox_layout = Layout(display='flex', flex_flow='column', align_items='center')

ui = widgets.VBox([html_label, Fo_slider_step3], layout=vbox_layout)
out = widgets.interactive_output(update_plots, {'selected_frequencies_part3': Fo_slider_step3})

# Display the slider
clear_output(wait=True)  # Clear the previous plot
display(ui,out)
<<<<<<< HEAD <<<<<<< HEAD
=======
>>>>>>> d546939e94d695c6e9c44c085bc5ec2e2e14440e ======= >>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

Part 4: Application B#

Write a Python spectral analysis function, similar to signal.welch(): it will accept as input a real signal vector as well as the sampling frequency, \(F_s\), and will plot the one-sided spectral density of the signal in the range \([0-F_s/2)\). The signal will be segmented into sections of length equal to the power of \(2\) closest to \(1/8\) of its total length, but not less than 256. The sections will overlap by \(50%\). The last section, if shorter than the others, will be ignored. The spectrum of each section will be calculated with FFT and the mean value of all sections will be taken. The function should be tested with the signal from example 1.1 and the result compared with that of signal.welch().

Code
# ==============================================================================
# Custom nextpow2 and pwelch functions
# ==============================================================================
def nextpow2(i):
        n = 1
        while n < i:
            n *= 2
        return n

def pwelch(x,Fs):                    
    Ts=1/Fs                    
    L=np.size(x)+1                 
    T=L*Ts                     
    N = 2^nextpow2(L)
    Fo=Fs/N                   
    f=np.arange(0,N)*Fo       
     
    window_size = nextpow2(np.size(x)/8)
    if (window_size<256):
        window_size=256
    windows = np.size(x)//(window_size//2)-1
    indexer = np.arange(window_size)[None, :] + (window_size//2)*np.arange(windows)[:, None]
    windowed_x = x[indexer]

    avg_pwr=0
    for window in windowed_x:
        window = window * np.hanning(np.size(window))
        L=np.size(window)+1                 
        T=L*Ts                     
        N = 2^nextpow2(L)
        Fo=Fs/N                   
        f=np.arange(0,N)*Fo
        window_fft=np.fft.fft(window,N)
        power=np.multiply(window_fft,np.conj(window_fft))/N/L
        avg_pwr=avg_pwr+power
    avg_pwr=avg_pwr/windows

    fig, ax = plt.subplots()
    fig.set_size_inches(18.5, 10.5)
    ax.plot(f[np.arange(0,N//2)],avg_pwr[np.arange(0,N//2)])
    ax.set(xlabel='Frequency (Hz)', ylabel='Power',
           title='Periodogram pwelch()')
    ax.grid()
    plt.show()
    
    return f[np.arange(0,N//2)], avg_pwr[np.arange(0,N//2)]

# ==============================================================================
# Test the functions & plot the results
# ==============================================================================   
Fs=______
f1,Pxx1 = pwelch(x,Fs)
f2,Pxx2 = signal.welch(x,fs=Fs)

fig, ax = plt.subplots()
fig.set_size_inches(18.5, 10.5)
ax.plot(f2,Pxx2)
ax.set(xlabel='Frequency (Hz)', ylabel='Power',
       title='Periodogram signal.welch()')
ax.grid()
plt.show()

% ==============================================================================
% Custom nextpow2 and pwelch functions
% ==============================================================================
function N = nextpow2(i)
    n = 1;
    while n < i
        n = n * 2;
    end
    N = log2(n); % This returns the exponent N such that 2^N = n, aligning with MATLAB's built-in nextpow2 behavior.
end

function [f, avgPwr] = pwelch(x, Fs)
    Ts = 1/Fs;
    L = length(x)+1;
    T = L * Ts;
    N = 2^nextpow2(L);
    Fo = Fs / N;
    f = (0:N-1) * Fo;
    
    windowSize = nextpow2(length(x)/8);
    if (windowSize < 256)
        windowSize = 256;
    end
    windows = floor(length(x) / (windowSize / 2)) - 1;
    indexer = bsxfun(@plus, (1:windowSize)', (windowSize / 2) * (0:windows-1));
    windowedX = x(indexer);
    
    avgPwr = 0;
    for i = 1:size(windowedX, 2)
        window = windowedX(:, i) .* hanning(size(windowedX, 1));
        L = length(window)+1;
        T = L * Ts;
        N = 2^nextpow2(L);
        Fo = Fs / N;
        f = (0:N-1) * Fo;
        windowFFT = fft(window, N);
        power = (windowFFT .* conj(windowFFT)) / (N * L);
        avgPwr = avgPwr + power;
    end
    avgPwr = avgPwr / windows;
    
    figure;
    set(gcf, 'Position', [100, 100, 1850, 1050]);
    plot(f(1:floor(N/2)), avgPwr(1:floor(N/2)));
    xlabel('Frequency (Hz)');
    ylabel('Power');
    title('Periodogram pwelch()');
    grid on;
    
    f = f(1:floor(N/2));
    avgPwr = avgPwr(1:floor(N/2));
end

% ==============================================================================
% Test the functions & plot the results
% ==============================================================================
Fs = _______; % Define your sampling frequency
x = randn(1, 10000); % Example data, replace with your actual data

[f1, Pxx1] = pwelch(x, Fs);
[f2, Pxx2] = pwelch(x, Fs, [], [], Fs); % MATLAB's built-in function

figure;
set(gcf, 'Position', [100, 100, 1850, 1050]);
plot(f2, Pxx2);
xlabel('Frequency (Hz)');
ylabel('Power');
title('Periodogram signal.welch()');
grid on;

def nextpow2(i):
        n = 1
        while n < i:
            n *= 2
        return n

def pwelch(x,Fs):                    
    Ts=1/Fs                    
    L=np.size(x)+1                 
    T=L*Ts                     
    N = 2^nextpow2(L)
    Fo=Fs/N                   
    f=np.arange(0,N)*Fo       
     
    window_size = nextpow2(np.size(x)/8)
    if (window_size<256):
        window_size=256
    windows = np.size(x)//(window_size//2)-1
    indexer = np.arange(window_size)[None, :] + (window_size//2)*np.arange(windows)[:, None]
    windowed_x = x[indexer]

    avg_pwr=0
    for window in windowed_x:
        window = window * np.hanning(np.size(window))
        L=np.size(window)+1                 
        T=L*Ts                     
        N = 2^nextpow2(L)
        Fo=Fs/N                   
        f=np.arange(0,N)*Fo
        window_fft=np.fft.fft(window,N)
        power=np.multiply(window_fft,np.conj(window_fft))/N/L
        avg_pwr=avg_pwr+power
    avg_pwr=avg_pwr/windows

    
    
    return f[np.arange(0,N//2)], avg_pwr[np.arange(0,N//2)]

# Function to update plots based on slider value
def update_plots(Fs):
   
    T = 1 / Fs  # Update sampling period
    t1 = np.arange(0, L) * T  # Update time vector
    
    # Recompute signal x with new sampling frequency
    last_x = np.sin(2 * np.pi * 30 * t1) + 0.8 * np.sin(2 * np.pi * 80 * (t1 - 2)) + np.sin(2 * np.pi * 60 * t1)
    
    # Compute pwelch
    f1, Pxx1 = pwelch(last_x, Fs)
    
    # Compute signal.welch
    f2, Pxx2 = signal.welch(last_x, fs=Fs)
    
    # Plot
    fig, axs = plt.subplots(2, 1, figsize=(18.5, 20))
    
    # Plot custom pwelch
    axs[0].plot(f1, Pxx1)
    axs[0].set(xlabel='Frequency (Hz)', ylabel='Power', title='Periodogram pwelch()')
    axs[0].grid()
    
    # Plot signal.welch
    axs[1].plot(f2, Pxx2)
    axs[1].set(xlabel='Frequency (Hz)', ylabel='Power', title='Periodogram signal.welch()')
    axs[1].grid()
    
    plt.tight_layout()

# Create slider for Fs
Fs_slider = widgets.IntSlider(
    value=500,
    min=100,
    max=2000, 
    step=100,
    description='Sampling Frequency (Fs):',
    layout=Layout(width='90%'),
    style={'description_width': 'initial'}, 
    continuous_update=False
)


html_label = widgets.HTML(
    value="""
    
    """
)

# Display the sliders and output
vbox_layout = Layout(display='flex', flex_flow='column', align_items='center')

ui = widgets.VBox([html_label, Fs_slider], layout=vbox_layout)
out = widgets.interactive_output(update_plots, {'Fs': Fs_slider})

# Display the slider
clear_output(wait=True)  # Clear the previous plot
display(ui,out)
<<<<<<< HEAD <<<<<<< HEAD
=======
>>>>>>> d546939e94d695c6e9c44c085bc5ec2e2e14440e ======= >>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4
Contents